#ifndef GAUSSLEGENDRERK_HPP
#define GAUSSLEGENDRERK_HPP

#include <iostream>
#include <vector>
#include <math.h>
#include <fstream>
#include "norm.hpp"
#include "Function.hpp"

using namespace std;

class GaussLegendreRK : public TimeIntegrator
{
private:
    int n_steps, p;
    double T, mu;
    vector<double> value;
    vector<vector<double>> result;
    Norm norm;

public:
    GaussLegendreRK() {}
    GaussLegendreRK(int n_, int _p, double time, vector<double> Val) : n_steps(n_), p(_p), T(time), value(Val) {}
    ~GaussLegendreRK() {}
    void solver();
    void error();
};

void GaussLegendreRK::solver()
{
    auto epsilon = 1e-9;
    ofstream outfile;
    outfile.open("ouput.m");
    outfile << "points = [";
    double k = 1.0 * T / n_steps;
    vector<double> U = value;
    outfile << U << endl;
    result.push_back(U);

    if (p == 1)
    {
        for (int i = 0; i < n_steps; i++)
        {
            auto y1 = U;
            while (true)
            {
                auto temp = f(U + (0.5 * k) * y1, i * k + 0.5 * k);
                if (norm.inftynorm(y1 - temp) < epsilon)
                    break;
                else
                    y1 = temp;
            }
            U = U + k * y1;
            outfile << U << endl;
            result.push_back(U);
        }
    }
    else if (p == 2)
    {
        for (int i = 0; i < n_steps; i++)
        {
            auto y1 = U;
            auto y2 = U;
            while (true)
            {
                auto temp1 = f(U + (0.25 * k) * y1 + ((3 - 2 * sqrt(3)) * k / 12) * y2, i * k + ((3 - sqrt(3)) * k / 6));
                auto temp2 = f(U + ((3 + 2 * sqrt(3)) * k / 12) * y1 + (0.25 * k) * y2, i * k + ((3 + sqrt(3)) * k / 6));
                if (norm.inftynorm(y1 - temp1) + norm.inftynorm(y2 - temp2) < epsilon)
                    break;
                else
                {
                    y1 = temp1;
                    y2 = temp2;
                }
            }
            U = U + k * (0.5 * y1 + 0.5 * y2);
            outfile << U << endl;
            result.push_back(U);
        }
    }
    else if (p == 3)
    {
        for (int i = 0; i < n_steps; i++)
        {
            auto y1 = U;
            auto y2 = U;
            auto y3 = U;
            while (true)
            {
                auto temp1 = f(U + (5 * k / 36) * y1 + ((2.0 / 9 - sqrt(15) / 15) * k) * y2 + ((5.0 / 36 - sqrt(15) / 30) * k) * y3, i * k + 0.1 * (5 - sqrt(15)) * k);
                auto temp2 = f(U + ((5.0 / 36 + sqrt(15) / 24) * k) * y1 + (2 * k / 9) * y2 + ((5.0 / 36 - sqrt(15) / 24) * k) * y3, i * k + 0.5 * k);
                auto temp3 = f(U + ((5.0 / 36 + sqrt(15) / 30) * k) * y1 + ((2.0 / 9 + sqrt(15) / 15) * k) * y2 + (5 * k / 36) * y3, i * k + 0.1 * (5 + sqrt(15)) * k);
                if (norm.inftynorm(y1 - temp1) + norm.inftynorm(y2 - temp2) + norm.inftynorm(y3 - temp3) < epsilon)
                    break;
                y1 = temp1;
                y2 = temp2;
                y3 = temp3;
            }
            U = U + (5 * k / 18) * y1 + (4 * k / 9) * y2 + (5 * k / 18) * y3;
            outfile << U << endl;
            result.push_back(U);
        }
    }
    outfile << "];" << endl;
    outfile << "scatter(points(:,1), points(:,2),\".\");" << endl;
    outfile.close();
}

void GaussLegendreRK::error()
{
    // cout << result[n_steps] << endl;
    // cout << result[0] << endl;
    double error = norm.Inftynorm(result[n_steps] - result[0]);
    cout << error << endl;
}

TimeIntegrator* createGaussLegend(int n_, int p_, double T_, vector<double> value_) { return new GaussLegendreRK(n_, p_, T_, value_); }
#endif